2020-12-14作成
2024-01-26更新
このドキュメントは、解剖学・骨学実習における計測・解析編の補足資料です。本実習では、霊長類の頭骨標本を使用し、形態計測、CT画像の処理、幾何学的形態測定を行います。これらの演習を通じて、基本的な技術の習得と関連概念の理解を目指します。
CT画像のセグメンテーションに特化した発展編はこちら。
ノギス、MicroScribe、一部の3Dスキャナ、デジタルカメラは容易に持ち運びできるので、たとえば博物館を渡り歩いてたくさんのデータを収集するといったことができる。CTやMRIは大掛かりな装置を必要とするが、標本内部の構造を観察することができる。画像データは各種デジタルレポジトリでの整備・公開が進んでおり、世界中の研究者によって有効利用されている。
形態画像データを収録するデジタルレポジトリの例
デジタルカメラで標本を撮像する場合は、視差によって生じる計測誤差に注意する必要がある。これは、標本とレンズの距離を近づけすぎないようにすること、カメラとレンズのセッティングをサンプル間で統一することなどにより、ある程度低減することができるようだ(Mullin & Taylor 2002)。
本実習では、ノギスとMicroScribeの計測を体験する。また、μCT(Bruker社製Skyscan 1275)の撮像を見学する(X線ガラスバッチを着用していない人はCT装置の操作はできない。自分で操作する場合は、事前に所定の講習を受けて登録する必要がある)。
CTは、管球でX線を発生させ、被写体に照射し、反対側に設置された検出器で捉える。このとき、管電圧はX線の強度に、管電流はX線の量に影響する。X線は、被写体によって吸収され、その密度と厚みに応じて減退するので、検出器では被写体のX線透過率の投影分布が取得される。これをあらゆる方向から行い、それらをコンピュータ上で再構成することで、被写体の三次元的構造を計算することができる。以下に、再構成の概念図を示す。
観測された投影データから未知数Xを推定する。上の単純な例では、線形連立方程式の解を求めることで未知数を特定できる。
\[ X_1+X_2=2\\ X_3+X_4=1\\ X_1+X_3=2\\ X_2+X_4=1\\ X_1+X_4=3 \]
しかし実際には、未知数は巨大になり、また測定値は誤差を含むので、この通りにはいかない。実際のCT装置では、計算負荷の低い、逆投影法を基本とする解析的方法が普及した。以下に、単純逆投影法の概略を示す。
単純逆投影法は、図に示されるように、被写体が存在しないところにも値が入るので、ボケが生じてしまい実用に適さない。このボケを低減させるために、実際のCT装置では、逆投影の前の投影データに(エッジを強調する)フィルタ補正を施す方法(フィルタ補正逆投影法)が一般的に取られてきた。当然、フィルタ(再構成関数)の種類によって再構成画像の印象は大きく変わり、エッジ効果の高いフィルタを適用するほど、空間分解能は高くなるがノイズが増える傾向にある(平尾 2008)。被写体の種類などに応じて、適切な再構成関数を選択する必要がある。
ただし、フィルタ補正を適切に施しても、アーチファクトを完全に消すことはできない。近年では、コンピュータの性能の向上に伴って、代数的・統計的再構成方法が採用されるようになり、ノイズやアーチファクトの低減といった、従来の解析的方法の欠点の解決が期待されている。
再構成されたデータは、通常、TIFFなどの一般的な画像形式もしくはDICOMという医用画像の標準規格をみたす形式で、2D画像のスタックという形で出力される。各ピクセルはX線吸収係数を換算したCT値を持ち、通常、水が0、空気が-1000、となるようにキャリブレーションされている(Hounsfield unit: HU)。CT値が大きいほど、密度が高く、X線吸収係数が大きいことを表す。
以下に示すように、さまざまなソフトウェアがある。
本実習では3D Slicer 5.6.0を使用する。
Digital Morphology Museum, KUPRIからダウンロードした以下のデータを使用する。2023年12月8日現在、セキュリティ上の問題により、外部からのアクセスは閉鎖されているため、データは授業で配布する。
| PRICT | ID | 学名 | 略称 | 和名 |
|---|---|---|---|---|
| 504 | WPK-Mimi | Pan troglodytes | Pt | チンパンジー |
| 949 | PRI-673 | Cebus albifrons | Ca | シロガオオマキザル |
| 1055 | PRI-4536 | Presbytis melalophos | Pm | クロカンムリリーフモンキー |
| 1095 | PRI-6093 | Varecia variegata | Vv | シロクロエリマキキツネザル |
| 1267 | PRI-6494 | Macaca fuscata | Mf | ニホンザル |
| 1287 | PRI-2516 | Hylobates lar | Hl | シロテテナガザル |
dicomの読み込み
2D表示
3D表示(ボリュームレンダリング)
3Dサーフェスモデルを作成するためには、まずセグメンテーション(画像の2値化)を行う必要がある。
本実習では、もっとも簡単なグローバル閾値法によるセグメンテーションを行う。今回、閾値は、見た目で判断して適当に設定するが、厳密な精度が要求される場合はhalf-maximum height (HMH)法(Spoor et al. 1993; Coleman and Colbert 2007)により求めるのが望ましい。たとえば、Coleman and Colbert (2007)では、ランダムに選択した10箇所においてHMHを算出し、それらの平均値をグローバル閾値としている。
HMHは、境界領域におけるCT値の最大値と最小値の平均値である。以下に例を示す。
空気と骨の境界領域に引いた線分上で、CT値の分布を示す。
ここで、最大値は1755、最小値は-1025なので、HMHはそれらの平均値として365となる。
HMH(365)をグローバル閾値をとしてセグメンテーションすると、以下の図のようになる。
仮に、閾値を-900とすると、以下の図のようになる。HMHのそれよりも骨領域が膨らんでいることがわかる。
一方、閾値を900とすると、以下の図のように、骨領域が細くなる。
このように、閾値は、高すぎると実際よりも痩せてしまったり穴が空いてしまったりし、逆に低すぎると実際よりも膨らんでしまったりアーチファクトが映り込んでしまったりする。このような性質を持つため、CT値のバラツキが大きい領域をセグメンテーションしようとする場合、グローバル閾値法は適さない。グローバル閾値法のほか、以下に示すような様々なセグメンテーションの方法が開発されており、データの性質や求められる精度に応じて使い分けるとよい。各手法の適用例や精度評価については、Engelkes (2020)、Rathnayaka et al. (2011)、Scherf & Tilgner (2009)、Ito (2019)などに解説がある。
グローバル閾値法によるセグメンテーション
3Dモデルの保存
ランドマークデータを取得するソフトウェアには、例えば以下のものがある。セミランドマークについては後述する。
本実習では再び3D Slicer 4.11を使用し、以下のランドマークのデータを取得する。
| Order | Landmark | Type |
|---|---|---|
| 1 | Prosthion | 2 |
| 2 | Glabella | 3 |
| 3 | Inion | 2 |
| 4 | Basion | 2 |
| 5 | Zygomaxillare (left) | 2 |
| 6 | Zygomaxillare (right) | 2 |
Bookstein (1991)は、ランドマークを3つのタイプに分類している。
- Type 1: discrete juxtapositions of tissues
- Type 2: maxima of curvature or other local morphogenetic processes
- Type 3: extremal points
幾何学的形態測定(次節で解説)では、ランドマークが相同であることを前提とする。したがって、確実に相同性を担保できるタイプ1のランドマークがもっとも望ましく、そうではないタイプ3のランドマークはあまり望ましくない。タイプ2のランドマークはそれらの中間的性質をもつ。
ランドマーキング
ランドマークデータの保存
fscvファイルの中を見てみよう。
head landmark/Pan_troglodytes.fcsv
## # Markups fiducial file version = 4.11
## # CoordinateSystem = LPS
## # columns = id,x,y,z,ow,ox,oy,oz,vis,sel,lock,label,desc,associatedNodeID
## 1,-2.0602643489837646,-87.3006591796875,-1103.945068359375,0,0,0,1,1,1,0,F-1,,
## 2,-0.43642905354499817,-30.19694709777832,-1016.7828979492188,0,0,0,1,1,1,0,F-2,,
## 3,2.1217479705810547,100.34082794189453,-1020.6141967773438,0,0,0,1,1,1,0,F-3,,
## 4,3.224494457244873,62.182159423828125,-1080.2637939453125,0,0,0,1,1,1,0,F-4,,
## 5,46.13127517700195,-23.04610824584961,-1078.1734619140625,0,0,0,1,1,1,0,F-5,,
## 6,-42.983951568603516,-23.148298263549805,-1081.05810546875,0,0,0,1,1,1,0,F-6,,
はじめの3行はコメントで、4行目から取得した順番に座標データが並ぶ。1列目は番号で、2〜3列目にXYZ座標データがミリ単位で記録されている。
次節では、このデータをRに読み込んで、幾何学的形態測定による解析を行う。
幾何学的形態測定(geometric morphometrics)では、古典的な形態測定法のように距離や角度に落とし込むことなく、座標データそのものを扱う。幾何学的形態測定では、拡大・縮小(scaling)、 回転(rotation)、移動(translation/centering)といった座標変換を施すことで、サンプル間の大きさ、方向、位置の違いを調整する。それにはいくつかの方法があるが、generalized Procrustes analysis (GPA)が一般的に使用される。このとき、大きさ(サイズ)は重心サイズ(Centroid size;重心から各標識点までの二乗距離の総和の平方根)で表される。整列後のデータは、Procrustes shapeと呼ばれ、これをシェイプ(形状)として利用する。
ここで注意が必要となるのは、整列後のデータの個別のランドマークに生物学的な意味を見出すのは困難ということである。整列後のデータはセットとして扱い、個別のランドマークを解析に供することは多くの場合ミスリーディングである(Cardini 2022)。
シェイプはサイズの成分を含まない(重心サイズが1に揃えられている)。ただしこれは、シェイプがサイズと無相関である(アロメトリー傾向がない)ということではない。
幾何学的形態測定については、Bookstein (1991)、Zelditch et al. (2012)、Mitteroecker & Gunz (2009)、Slice (2007)、Adams et al. (2004)などの参考文献がある。
今回の実習では、簡単のため6サンプルのみを解析に用いるが、ランドマーク数に対してサンプルサイズが少なすぎると統計的な問題となる場合があるので、本来は望ましくない。実際の研究では、十分なサンプルサイズを確保するようにしたい。
手始めに、Pan troglodytesのランドマークデータを読み込んでみる。
library(tidyverse)
read_csv("landmark/Pan_troglodytes.fcsv",
skip = 3, # 最初の3行のコメント行をスキップ
col_names = FALSE)[,2:4] # 座標データが含まれる2,3,4列を抽出
## # A tibble: 6 × 3
## X2 X3 X4
## <dbl> <dbl> <dbl>
## 1 -2.06 -87.3 -1104.
## 2 -0.436 -30.2 -1017.
## 3 2.12 100. -1021.
## 4 3.22 62.2 -1080.
## 5 46.1 -23.0 -1078.
## 6 -43.0 -23.1 -1081.
次に、関数を作成し、指定したフォルダに含まれるすべてのデータを一挙に読み込み配列(array)に格納する。
read_fcsv <- function(path, k = 3, p = 6){
list_data <- lapply(
list.files(path, full.names = TRUE),
function(x){read_csv(x, skip = 3, col_names = FALSE, col_types = cols())}[,2:4]%>%as.matrix)
# 配列に格納する
array_data <- array(as.numeric(unlist(list_data)), dim=c(p, k, length(list_data)))
# 名前をつける
names <- tools::file_path_sans_ext(list.files(path, pattern = "*.csv"))
dimnames(array_data)[[3]] <- names
return (array_data)
}
raw_data <- read_fcsv("landmark")
dimnames(raw_data)[[3]] <- c("Ca", "Hl", "Mf", "Pt", "Pm", "Vv") # 種名が長いので略称に
raw_data
## , , Ca
##
## [,1] [,2] [,3]
## [1,] -0.1149868 -35.12373 -1054.850
## [2,] -0.6637123 -17.67938 -1023.330
## [3,] -0.1195280 53.23666 -1045.329
## [4,] -0.1068767 23.19316 -1060.824
## [5,] 19.9068317 -13.63381 -1053.072
## [6,] -20.5708523 -13.87484 -1053.465
##
## , , Hl
##
## [,1] [,2] [,3]
## [1,] 1.4124396 49.33804 -1338.383
## [2,] -0.2931501 30.91640 -1305.041
## [3,] -0.8856894 -54.35939 -1319.973
## [4,] 2.1998587 -23.15752 -1341.715
## [5,] -23.5790615 21.08697 -1340.217
## [6,] 23.8187637 21.20558 -1338.248
##
## , , Mf
##
## [,1] [,2] [,3]
## [1,] -0.8745944 -51.14996 -1061.536
## [2,] -0.9600250 -17.43713 -1020.646
## [3,] -1.9526790 65.18864 -1053.153
## [4,] -1.5969342 28.84805 -1069.347
## [5,] 27.5209427 -15.30799 -1064.343
## [6,] -27.5567760 -15.59646 -1066.165
##
## , , Pt
##
## [,1] [,2] [,3]
## [1,] -2.0602643 -87.30066 -1103.945
## [2,] -0.4364291 -30.19695 -1016.783
## [3,] 2.1217480 100.34083 -1020.614
## [4,] 3.2244945 62.18216 -1080.264
## [5,] 46.1312752 -23.04611 -1078.173
## [6,] -42.9839516 -23.14830 -1081.058
##
## , , Pm
##
## [,1] [,2] [,3]
## [1,] -2.2882662 -38.44722 -1014.9694
## [2,] -0.9373286 -16.93668 -985.1452
## [3,] -0.3596495 54.95614 -1010.1675
## [4,] -0.9612941 25.87814 -1026.7108
## [5,] 24.1622162 -19.22473 -1017.1849
## [6,] -27.1201286 -17.04971 -1016.9243
##
## , , Vv
##
## [,1] [,2] [,3]
## [1,] 0.1620082 -52.117371 -1306.246
## [2,] -0.9636781 -15.641431 -1282.679
## [3,] -2.3839674 53.863892 -1289.016
## [4,] -0.2980522 41.737179 -1307.885
## [5,] 20.6963501 -3.194123 -1302.389
## [6,] -21.0222626 -3.935471 -1304.899
読み込んだ座標データを表示してみよう。
library(geomorph)
plotAllSpecimens(raw_data)
CT撮像時に標本の位置を揃えていないため、とくにZ方向の位置がバラバラである。方向は概ね揃えてあるが、厳密には揃っていない。当然、大きさもバラバラである。小さな黒い点は平均形状を示す。
そこで、Procrustes整列を行う。
gpa_data <- gpagen(raw_data, print.progress = FALSE)
plotAllSpecimens(gpa_data$coords)
サンプル間の大きさ、方向、位置の違いが調整され、整列されていることがわかる。
計測には誤差がつきものである。主解析の前に、誤差の程度を評価することが重要である。
同じ標本を2度ずつ計測したとする。2回目の計測データを読み込む。
raw_data2 <- read_fcsv("landmark2")
dimnames(raw_data2)[[3]] <- c("Ca", "Hl", "Mf", "Pt", "Pm", "Vv")
1回目の計測データと比較してみよう。
plotRefToTarget(raw_data[,,1], raw_data2[,,1], method = "point") # 一つ目の標本
次に、Procrustes ANOVA(分散分析)という方法で、誤差の程度を定量的に評価する。
combind_data <- arrayspecs(
A = rbind(
two.d.array(raw_data),
two.d.array(raw_data2)
),
p = 6,
k = 3
)
cranium_sym <- bilat.symmetry(
A = combind_data,
ind = dimnames(combind_data)[[3]],
replicate = c(rep(1, 6), rep(2, 6)),
object.sym = TRUE,
land.pairs = matrix(c(5, 6), nrow = 1, ncol = 2),
print.progress = FALSE
)
summary(cranium_sym)
##
## Call:
## bilat.symmetry(A = combind_data, ind = dimnames(combind_data)[[3]],
## replicate = c(rep(1, 6), rep(2, 6)), object.sym = TRUE, land.pairs = matrix(c(5,
## 6), nrow = 1, ncol = 2), print.progress = FALSE)
##
##
## Symmetry (data) type: Object
##
## Type I (Sequential) Sums of Squares and Cross-products
## Randomized Residual Permutation Procedure Used
## 1000 Permutations
##
## Shape ANOVA
## Df SS MS Rsq F Z Pr(>F)
## ind 5 0.32999 0.065998 0.98406 162.6990 1.96297 0.023 *
## side 1 0.00070 0.000698 0.00208 1.7199 0.72716 0.245
## ind:side 5 0.00203 0.000406 0.00605 1.8592 2.09381 0.021 *
## ind:side:replicate 12 0.00262 0.000218 0.00781
## Total 23 0.33533
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
今回のデータでは、個体差(種間差)に対して、計測誤差は、極めて小さいことがわかる。しかし、種内変異や個体内変異(左右差)などの微小な変異を扱うような研究では、計測誤差は無視できない程度になる可能性があるため、注意が必要である。
2回の計測の平均値を取ることで、計測誤差による偏りを軽減することができる。左右対称標本の場合、しばしば、左右対称形状成分と非左右対称形状成分に分けて解析することがある。左右差を興味の対象としない場合は、左右対称形状成分を解析に供すれば良いだろう。
cranium_sym$symm.shape
## , , Ca
##
## [,1] [,2] [,3]
## [1,] 0.4136886 -2.399107e-09 -0.09055265
## [2,] 0.2065382 -3.252927e-08 0.30117827
## [3,] -0.6451545 1.133182e-08 0.05773511
## [4,] -0.2915072 2.017206e-08 -0.14069160
## [5,] 0.1582174 2.376169e-01 -0.06383454
## [6,] 0.1582174 -2.376169e-01 -0.06383459
##
## , , Hl
##
## [,1] [,2] [,3]
## [1,] 0.4327259 -6.667695e-09 -0.04789490
## [2,] 0.2176267 -1.341479e-08 0.27931840
## [3,] -0.6351050 9.995024e-09 0.06478416
## [4,] -0.3012089 1.025778e-08 -0.13591751
## [5,] 0.1429807 2.442357e-01 -0.08014507
## [6,] 0.1429807 -2.442357e-01 -0.08014508
##
## , , Mf
##
## [,1] [,2] [,3]
## [1,] 0.4649953 4.887770e-09 -0.08221865
## [2,] 0.1669177 -2.718378e-08 0.32290939
## [3,] -0.6122245 -2.252171e-09 0.05791856
## [4,] -0.2820839 1.025574e-08 -0.11113496
## [5,] 0.1311977 2.535221e-01 -0.09373715
## [6,] 0.1311977 -2.535221e-01 -0.09373719
##
## , , Pt
##
## [,1] [,2] [,3]
## [1,] 0.5130677 1.687515e-08 -0.10376616
## [2,] 0.1041127 -1.705663e-08 0.28535792
## [3,] -0.5890496 -1.812760e-08 0.10092188
## [4,] -0.3122893 4.938586e-09 -0.16446396
## [5,] 0.1420793 2.405358e-01 -0.05902482
## [6,] 0.1420793 -2.405358e-01 -0.05902485
##
## , , Pm
##
## [,1] [,2] [,3]
## [1,] 0.4001911 -7.403814e-09 -0.06704527
## [2,] 0.1899733 -1.802722e-08 0.28262999
## [3,] -0.6175216 1.322295e-08 0.06559211
## [4,] -0.3164259 1.463837e-08 -0.13750531
## [5,] 0.1718915 2.782285e-01 -0.07183575
## [6,] 0.1718915 -2.782285e-01 -0.07183578
##
## , , Vv
##
## [,1] [,2] [,3]
## [1,] 0.60120675 2.928514e-08 -0.02683340
## [2,] 0.18687293 9.577850e-09 0.19636092
## [3,] -0.55321784 -2.649713e-08 0.06119630
## [4,] -0.39661730 -1.978524e-08 -0.13748963
## [5,] 0.08087772 2.282363e-01 -0.04661709
## [6,] 0.08087774 -2.282363e-01 -0.04661709
今回使用した標本は欠けているところがなかったので欠損値はないが、実際には一部が欠けている標本は珍しくない。また、化石標本に至っては完全なものの方が少ないだろう。しかし、幾何学的形態測定は欠損データを扱うことができない。したがって、欠損ランドマークを含む標本は予め解析から除外するか、もしくは欠損ランドマークを補完するする必要がある。
仮に、標本#1のランドマーク#1の座標データを取得することができなかったとする。ここでは、薄板スプライン関数(TPS)を用いた補完を行う。TPSにより、欠損のない参照データを欠損を含むデータに重ね合わせることで、欠損ランドマークの座標を推定する。
raw_data_missing <- raw_data
raw_data_missing[1,,1] <- NA #標本#1のランドマーク#1を欠損データとする
raw_data_tps <- estimate.missing(raw_data_missing, method = "TPS")
plotRefToTarget(raw_data[,,1], raw_data_tps[,,1], method = "point")
今回用いたデータセットの場合、実際のランドマークと補完されたランドマークは大きく異なっていることが分かる。欠損ランドマークを持つ個体に似た個体がデータセットに含まれない場合は注意が必要である。
TPSの他にも、回帰、ベイズ主成分分析、平均値置換、などによる欠損値補完の方法がある。それぞれの利点と欠点、注意点については、Arbour et al. (2014)に詳しい。
整列後のデータに対して、主成分分析を実行し、第一主成分と第三主成分の散布図を描く。
pca_result <- gm.prcomp(cranium_sym$symm.shape)
library(ggplot2)
pca_result$x %>%
as_tibble() %>%
mutate(Species = dimnames(raw_data)[[3]]) %>%
ggplot(aes(x = Comp1, y = Comp2, label = Species)) +
geom_point() +
geom_text(hjust = 0, nudge_x = 0.005)
各主成分の寄与率(proportion of variance)を表示する。
tibble(
PC = 1:5,
Proportion_of_variance = pca_result$sdev^2/sum(pca_result$sdev^2),
Cummurative_proportion = cumsum(pca_result$sdev^2/sum(pca_result$sdev^2))
) %>%
mutate(across(2:3, ~ round(., 2)))
## # A tibble: 5 × 3
## PC Proportion_of_variance Cummurative_proportion
## <int> <dbl> <dbl>
## 1 1 0.71 0.71
## 2 2 0.17 0.89
## 3 3 0.07 0.95
## 4 4 0.04 0.99
## 5 5 0.01 1
第一主成分と第二主成分で、全主成分のうち88%の分散を説明している。
第一主成分の形状変化を可視化する。平均形状から第一主成分の最大値への変化を示す。
plotRefToTarget(pca_result$shapes$shapes.comp1$min, pca_result$shapes$shapes.comp1$max, method = "vector", label = TRUE)
第一主成分が大きいほど、全体的に細長い傾向にあるようだ。
続いて、第二主成分。
plotRefToTarget(pca_result$shapes$shapes.comp2$min, pca_result$shapes$shapes.comp2$max, method = "vector", label = TRUE)
第二主成分が大きいほど、相対的に顔面高が大きい傾向にあるようだ。
形状データの距離行列を用いて、系統関係を推定してみよう。
手始めに、CaとHlの間のユークリッド距離を計算する。 \[d = \sqrt{\sum(X_i - Y_i)^2} \]
sqrt(
sum(
(pca_result$x[1, ] -
pca_result$x[2, ]
)^2
)
)
## [1] 0.06430009
これは、dist関数を用いることで、一挙に計算できる。
dist(pca_result$x)
## Ca Hl Mf Pt Pm
## Hl 0.06430009
## Mf 0.10268928 0.09519065
## Pt 0.16548752 0.16668579 0.13136411
## Pm 0.08130329 0.08285905 0.11848060 0.17093575
## Vv 0.28783854 0.25090738 0.25714266 0.21739240 0.28792127
これで距離行列の準備ができた。では、まずは、もっとも単純な方法の一つである、UPGMA(Unweighted pair group method with arithmetic mean)法で系統を推定してみよう。
このうち、最小の距離をもつ組み合わせは、CaとHl(0.0643001)である。これらを一つにまとめ、クラスタとし、距離行列をアップデートする。ここで、各要素とクラスタとの距離は、そのクラスタに含まれる全ての要素との距離の平均とする。
## Mf Pt Pm Vv
## Pt 0.13136411
## Pm 0.11848060 0.17093575
## Vv 0.25714266 0.21739240 0.28792127
## (Ca,Hl) 0.09893996 0.16608666 0.08208117 0.26937296
枝長は距離の半分として0.03215となる。
アップデートされた距離行列の中で、最小の距離をもつ組み合わせは、Pmと(Ca,Hl)(0.0820812)である。これらをまとめると次のようになる。
## Mf Pt Vv
## Pt 0.1313641
## Vv 0.2571427 0.2173924
## (Pm,(Ca,Hl)) 0.1054535 0.1677030 0.2755557
以下同様に、最小の距離をもつものから、クラスタにまとめ上げていく。
## Pt Vv
## Vv 0.2173924
## (Mf,(Pm,(Ca,Hl))) 0.1586183 0.2709525
UPGMAは、hclust関数を使って計算することもできる。
library(ape)
plot(
as.phylo(
hclust(dist(pca_result$x), method = "average")
)
)
axisPhylo()
UPGMAは進化速度の一定性を仮定しており、有根系統樹が得られる。
一方、進化速度一定の仮定を要しない方法としては、近隣結合法(neighbor joining)がある。
nj_morph <- nj(dist(pca_result$x))
plot(nj_morph ,type="unrooted")
Vvを外群とした時の有根系統樹(頭蓋形状)
nj_morph_rooted <- root(nj_morph, "Vv")
plot(nj_morph_rooted)
axisPhylo()
UPGMAと概ね同じ分岐関係が得られるが、Hl、Ca、Pmの分岐関係は異なっている。
では、頭蓋形状のデータから推定した系統は、系統関係を正しく推定できているだろうか。 DNAの塩基配列に基づく系統推定の結果と比較してみよう。 10kTreesからアライメント済のデータ [Final dataset (NEXUS format)] をダウンロードして、Rに読み込み、近隣結合法により系統を推定する。
# 10kTreesのnexus dataの読み込み
nexus_10ktrees <- read.nexus.data("10kTrees/10kTrees_finalFile_version3.nex")
# 形態データと同じ6種のデータを抽出
nexus_sub <- nexus_10ktrees[c("Cebus_albifrons", "Macaca_fuscata", "Pan_troglodytes_verus", "Presbytis_melalophos", "Varecia_variegata_variegata", "Hylobates_lar")]
# DNAbin形式に変換
dnabin_sub <- nexus2DNAbin(nexus_sub)
# 名前を略称に
names(dnabin_sub) <- c("Ca", "Mf", "Pt", "Pm", "Vv", "Hl")
nj_dna <- nj(dist.dna(dnabin_sub))
plot(nj_dna,type="unrooted")
Vvを外群とした時の有根系統樹(DNA塩基配列)
nj_dna_rooted <- root(nj_dna, "Vv")
plot(nj_dna_rooted)
axisPhylo()
DNA塩基配列に基づく系統樹では、旧世界ザル類のMfとPm、ヒト上科のPtとHlがそれぞれクラスターを形成し、その外側に新世界ザル類のCaが来ており、枝長はともかく分岐関係は正しく推定できているように見える。一方、形状に基づく系統樹は、まったく間違った分岐関係になっている。(少なくとも今回のデータで表される)頭蓋形状の変異は、系統関係を反映しないようだ。
幾何学的形態測定のデータに基づく系統推定の方法論については、未だ活発に議論が続いている。よくある落とし穴は、一部の主成分の得点やProcrustes整列後の座標データを系統解析の形質として用いるというもので、統計的な問題があるようだ(Adams et al. 2011)。本実習で紹介した、Procrustes距離行列に基づくUPGMAや近隣結合法のほか、最節約法による推定方法も開発されている(Catalano et al. 2011)。
最後に、形状の主成分得点の散布図に、DNA塩基配列に基づく系統樹をプロットしてみよう。
pca_result <- gm.prcomp(gpa_data$coords, phy = nj_dna)
plot(pca_result, phylo = TRUE)
これはphylomorphospaceと呼ばれる解析手法で、主成分で表される形状空間において、系統発生の過程でどのような変化があったのかが可視化される。これは単に主成分空間に後から系統樹をプロットしたもので、主成分自体は先に示したものと同一である。
解剖学的な特徴点のない部位の形状はどうしたら評価できるだろうか。一つには、セミランドマークを使った方法が考えられる。以下に、その方法を示す。
正中面における脳頭蓋の輪郭の形状を評価するとしよう。
ランドマーキング
ランドマークと同様の方法で保存し、Rに読み込む。
raw_sl <- read_fcsv("semilandmark", k = 3, p = 15)
new_data <- arrayspecs(
cbind(
two.d.array(raw_data[-c(2,3),,]),
two.d.array(raw_sl)
),
k = 3,
p = 19
)
Procrustes整列をする際に、セミランドマークをスライドさせる。これは、輪郭の接線方向の変動(生物学的な意味を持たない)を除去するための処置である。この際、屈曲エネルギーを最小化させる方法と、Procrustes距離を最小化する方法があるが、今回は前者を利用する。それぞれの利点や欠点、注意点などは、Perez et al. (2006)に詳しい。
gpa_data2 <- gpagen(A = new_data, curves = define.sliders(1:15), print.progress = FALSE)
左右対称成分を主成分分析に供し、その結果を示す。
cranium_sym2 <- bilat.symmetry(A = gpa_data2$coords, ind = dimnames(gpa_data2$coords)[[3]], land.pairs = matrix(c(3,4),ncol = 2), object.sym = TRUE, print.progress = FALSE)
pca_result2 <- gm.prcomp(cranium_sym2$sym)
pca_result2$x %>%
as_tibble() %>%
mutate(Species = dimnames(raw_data)[[3]]) %>%
ggplot(aes(x = Comp1, y = Comp2, label = Species)) +
geom_point() +
geom_text(hjust = 0, nudge_x = 0.005)
plotRefToTarget(pca_result2$shapes$shapes.comp1$min, pca_result2$shapes$shapes.comp1$max, method = "vector")
plotRefToTarget(pca_result2$shapes$shapes.comp2$min, pca_result2$shapes$shapes.comp2$max, method = "vector")
今回、セミランドマークを用いて曲線の形状を評価したが、曲面の形状を評価することも、同様の方法で可能である(Gunz et al. 2005)。